2105.10508v2 [astro-ph.GA] 6 Aug 2021 


e 
. 


arXiv 


MNRAS 000, 1-20 (2021) Preprint 9 August 2021 Compiled using MNRAS IATEX style file v3.0 


A universal relation between the properties of supermassive black holes, 
galaxies, and dark matter halos 


A. Marasco,!* G. Cresci,! L. Posti? F. Fraternali,? F. Mannucci,! A. Marconi, ^ F. Belfiore, ! S. M. Fall? 
VINAF - Osservatorio Astrofisico di Arcetri, Largo E. Fermi 5, 50127, Firenze, Italy 

? Université de Strasbourg, CNRS UMR 7550, Observatoire astronomique de Strasbourg, 11 rue de l'Université, 67000 Strasbourg, France 

3 Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands 

^Dipartimento di Fisica e Astronomia, Università di Firenze, Via G. Sansone 1, 50019, Sesto Fiorentino (Firenze), Italy 

5Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA 


Accepted XXX. Received Y Y Y; in original form ZZZ 


ABSTRACT 

We study the relations between the mass of the central black hole (BH) Mgu, the dark matter halo mass My, and the stellar-to-halo 
mass fraction f, cc M, / My in a sample of 55 nearby galaxies with dynamically measured Mgy > 106 Mo and My > 5x10!! Mo. 
The main improvement with respect to previous studies is that we consider both early- and late-type systems for which Mp is 
determined either from globular cluster dynamics or from spatially resolved rotation curves. Independently of their structural 
properties, galaxies in our sample build a well defined sequence in the Mpu-Mnp- f. space. We find that: (1) Mp and Mpu strongly 
correlate with each other and anti-correlate with fx; (ii) there is a break in the slope of the Mpu-Miy relation at Mp of 10? Mo, 
and in the f,-Mgy relation at Mgy of ~ 107—108 Mo; (iii) at a fixed Mpy, galaxies with a larger f, tend to occupy lighter halos 
and to have later morphological types. We show that the observed trends can be reproduced by a simple equilibrium model in 
the ACDM framework where galaxies smoothly accrete dark and baryonic matter at a cosmological rate, having their stellar and 
black hole build-up regulated both by the cooling of the available gas reservoir and by the negative feedback from star formation 
and active galactic nuclei (AGN). Feature (ii) arises as the BH population transits from a rapidly accreting phase to a more gentle 
and self-regulated growth, while scatter in the AGN feedback efficiency can account for feature (iii). 
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1 INTRODUCTION feedback depends on the BH accretion rate Mgy. Thus, in this pic- 
ture, galaxies can be fully described by the co-evolution among their 
stellar, BH and dark matter contents, whose growths are intertwined. 
A direct consequence of this evolutionary scenario is that, at any 
redshift, Mp, M, and Mgg are expected to be related to each other. 
A careful characterisation of the relations between these three quan- 
tities in nearby systems can give fundamental clues on the parameters 
that regulate galaxy evolution, and is the subject of this work. 


In a simplified theoretical framework, the build-up of galaxies can 
be thought as resulting from the competition between the ‘positive’ 
process of cooling and gravitational collapse of gas within the po- 
tential wells provided by dark matter halos (White & Rees 1978), 
and a series of ‘negative’ mechanisms such as gas heating and sub- 
sequent expulsions caused by feedback from star formation (Larson 
1974; Dekel & Silk 1986) and active galactic nuclei (AGN; Silk & 
Rees 1998; Harrison 2017). This competition between inflows and 
outflows or, alternatively, cooling and heating, ultimately regulates 
the galaxy gas reservoir out of which stars form and super-massive 
black holes (BHs) grow. This simple framework is at the basis of sev- 
eral successful theoretical models of galaxy formation and evolution 
(e.g. Somerville et al. 2008; Bouché et al. 2010; Lilly et al. 2013; 
Behroozi et al. 2019). 

The picture described above suggests the existence of three ‘lead- 
ing characters' playing a major role in the evolution of a galaxy, 
namely its dark matter halo, its stellar component and its BH. Their 
masses (My, My and Mpy) and growth rates are closely related to 
the positive and negative processes discussed: gas accretes onto ha- 
los at rates of fp My (fo being the Universal baryon fraction), stellar 
feedback depends on the star formation rate Mx (or SFR) and AGN 


Most studies in the literature have focused on the correlation be- 
tween pairs of this leading trio. In particular, the relation between 
My and M, or, equivalently, between Mp and the galaxy global star 
formation efficiency f, = M,/ fy My, has received a lot of atten- 
tion from both the observational and the theoretical communities. 
From the theoretical side, this so-called stellar-to-halo mass relation 
(SHMR; for a recent review see Wechsler & Tinker 2018) is com- 
monly probed via a semi-empirical technique known as abundance 
matching, which relates galaxies to halos by matching the observed 
stellar mass function to the theoretical halo mass function obtained 
from cosmological N-body simulations, assuming that stellar mass 
increases monotonically with the mass of the host halo. Different 
abundance matching studies (e.g. Vale & Ostriker 2004; Behroozi 
et al. 2010; Moster et al. 2013; Kravtsov et al. 2018) all point towards 
a scenario where f, is maximal in galaxies with Mp ~ 101? Mo (or 
My ~ 5x 101 Mo) and rapidly decreases at lower and higher Mp, 
* E-mail: antonino.marasco G inaf.it which is traditionally interpreted as evidence for negative feedback 
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from star formation in the low-mass regime and from AGN activity 
in the high-mass one. Observationally, the SHMR can be probed via 
different methods such as galaxy-galaxy weak lensing (Mandelbaum 
et al. 2006; Leauthaud et al. 2012), satellite or globular cluster kine- 
matics (More et al. 2011; van den Bosch et al. 2019), internal galaxy 
dynamics (Cappellari et al. 2013; Read et al. 2017) or a combination 
of these (Dutton et al. 2010). Recent studies from Posti & Fall (2021) 
and Posti et al. (20192) have found strong evidence for a difference 
in the SHMR of early- and late-type systems, with the former fol- 
lowing the standard predictions from abundance matching models 
while the latter showing a monotonically increasing fx as a function 
of mass, with no distinctive ‘peak’ in f,. This discrepancy results 
in a substantially different f, associated to the two galaxy types at 
My ~ 10!! Mo (that is, where most massive spirals are observed), 
suggesting the existence of different pathways for the stellar mass 
build-up in early- and late-type systems. 

The existence of empirical correlations between the mass of the 
supermassive BH and the properties of the host galaxy bulge (like 
its mass Mpylge and velocity dispersion c") has been largely explored 
in the literature (for an in-depth review see Graham 2016) and is 
now well established (e.g. Magorrian et al. 1998; Marconi & Hunt 
2003; Kormendy & Ho 2013; Saglia et al. 2016; de Nicola et al. 
2019). The Mgy-o relation, in particular, is amongst the tightest 
ones, with a vertical intrinsic scatter of ~ 0.3 dex, which is often 
interpreted as evidence for co-evolution between the BH and the host 
bulge resulting from AGN feedback (King 2003; Somerville et al. 
2008; King & Pounds 2015). Mergers can also play a role (Peng 
2007; Jahnke & Macció 2011), and there is evidence that the slope 
and normalisation of the relation between Mgy and Mpulge (or 0) 
depend on the galaxy structural parameters such as the bulge-to-total 
ratio and Sersic index (Graham & Scott 2013; Scott et al. 2013; Sahu 
et al. 2019a,b). Spiral galaxies also exhibit a very tight (intrinsic 
scatter of ~ 0.3 dex) correlation between Mpy and the spiral arm 
pitch angle (Seigar et al. 2008; Berrier et al. 2013; Davis et al. 2017), 
which is also supported by cosmological hydrodynamical simulations 
(Mutlu-Pakdil et al. 2018). 

Surprisingly, however, the relation between Mpg and the total 
stellar mass of the host M, has been explored much less in the 
literature and with somewhat contradictory results, ranging from 
being non-existent (Kormendy & Gebhardt 2001) to being as tight 
as the Mpy-Mbuige relation (Lasker et al. 2014). McConnell & Ma 
(2013) and Davis et al. (2018) gave more moderate views on the 
subject, showing the existence of a positive correlation between Mgy 
and M, as expected from a co-evolution scenario (e.g. Bower et al. 
2017, hereafter B17), although with a larger scatter with respect to 
the relations with the bulge properties. 

Several theoretical models have suggested that the dynamically- 
dominant component of a galaxy, i.e. its dark matter, should dictate 
the formation of BHs (e.g. Loeb & Rasio 1994; Haehnelt et al. 1998; 
Booth & Schaye 2010, B17) as a direct consequence of the physics 
of gas accretion onto halos. On the observational side, the pioneer- 
ing works by Whitmore et al. (1979) and Whitmore & Kirshner 
(1981) highlighted the existence of a relation between the galaxy 
rotational velocity vrot, traced by H1 kinematics, and the central ve- 
locity dispersion of the stellar component. However, it was Ferrarese 
(2002) who firstly interpreted this relation as evidence for a correla- 
tion between Mpy and Mj, since dark matter dominates the galaxy 
kinematics at large radii. The study of Ferrarese and later develop- 
ments (e.g. Pizzella et al. 2005; Volonteri et al. 2011) were criticised 
by Kormendy & Bender (2011), who showed that the correlation 
was apparent only in galaxies hosting a classical bulge, blaming the 
‘rotation curve conspiracy’ (e.g. van Albada & Sancisi 1986) as a 
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possible culprit for the observed trend. Sabra et al. (2015) used a 
sample of 53 galaxies of different morphological types with direct 
(dynamical) measurements of Mpy and only found evidence for an 
extremely weak correlation between vrot and the BH mass. However, 
the sample of Sabra et al. (2015) was later expanded by Davis et al. 
(2019b) and Smith et al. (2021), who concluded that a MBH-Vrot 
relation for spiral galaxies exists, consistent with expectations from 
the joint Mpy — M, and Tully & Fisher (1977) relations. Similar 
results were derived by Robinson et al. (2021) using a sample of 
24 systems with Mgy measurements from reverberation mapping. 
One of the limitations of these studies is that they often make use 
of Wso, the line-width of the integrated velocity profile from Hı or 
CO emission-line measurements, as a proxy for Vrot, which may lead 
to spurious results in cases where the rotation curve declines in the 
inner regions or the gas is not sufficiently extended in radius (e.g. 
Brook et al. 2016; Ponomareva et al. 2017). 

These considerations indicate that, while several fragmented 
pieces of evidence for a co-evolution between stars, dark matter and 
BHs in galaxies exist, steps need to be taken in order to build a more 
coherent observational picture, preparatory for constraining our the- 
oretical understanding of galaxy evolution. The goal ofthis study is to 
provide important steps in this direction. On the one hand, we aim to 
clarify the relationship between Mpy, My and M, (or rather, we pre- 
fer to focus on fy instead of M.) from observational data. The main 
improvement with respect to previous works is that, from the large 
pool of systems with dynamical Mgg estimates, we select a suitable 
sub-sample with dynamical measurements of Mp coming, for the 
vast majority of objects, from either globular cluster kinematics (for 
galaxies of earlier Hubble types) or spatially resolved rotation curves 
from interferometric Hı data (for galaxies of later Hubble types). 
Unlike studies that use Ha or CO data, which are limited to the inner 
regions of the galaxy, or H1 line-width data from single-dish tele- 
scopes, which lack spatial resolution, the measurements used in this 
work allow to trace galaxy dynamics up to very large distances from 
the galaxy centres (typically ~ 50 kpc) and are better suited to model 
the dark and luminous matter distribution in galaxies. In addition, 
we embed our observational results within a more general theoret- 
ical framework, showing that the observed relations are consistent 
with simple evolutionary models in ACDM where the stellar and BH 
mass built-up are regulated by the competition between positive and 
negative mechanisms discussed at the beginning of this Section. 

This paper is organised as follows. Our sample of nearby galax- 
les is described in Section 2 and the resulting scaling relations are 
presented in Section 3. Section 4 is dedicated to the building and 
the application of our model of galaxy evolution. The limitations of 
our model and the comparison with previous works are discussed in 
Section 5. Finally, Section 6 presents a summary and the conclusions 
of this study. 


2 GALAXY SAMPLE 


While several works have focused on the determination of galaxy 
BH, halo and stellar masses separately, we are not aware of a com- 
prehensive study where these three quantities have been derived 
simultaneously in a homogeneous, self-consistent fashion. The main 
reasons for this are the very different physical scales associated to 
these masses, whose measurements requires very diverse data sets. 
Therefore, the galaxy sample used in this work is derived from a 
combination of different datasets, using a variety of approaches. We 
show that, in spite of such diversity, a coherent picture emerges. All 
halo masses presented in this study are computed within the radius 
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Table 1. Main properties for the sample of 55 nearby galaxies studied in this work. 


Galaxy T-type  logjg(Magu/ Mo)  logj9(Mn/ Mo) log ig Cf) Ref. for Mau Mp method Ref. for My 

a) Q) E) (4) 6) (6) (7) (8) 

Milky Way 4.0 6.63 + 0.04 12.11 + 0.10 —0.68 + 0.18 KH13 GCs Posti & Helmi (2019) 

M31 3.0 8.15 +0.16 12.08 + 0.30 —0.24 + 0.32 KH13 RC (H1) Corbelli et al. (2010) 

M66 3.1 6.93 + 0.05 12.09 + 0.24 —0.62 + 0.28 S16 RC (Ho) Chemin et al. (2003) 

M81 2.4 7.81 x 0.13 12.01 + 0.24 —0.61 x 0.28 S16 RC (H1) de Blok et al. (2008) 

Centaurus A -2.1 7.75 + 0.08 12.34 + 0.24 —0.81 + 0.28 S16 RC (H1) vanGorkom et al. (1990) 

Circinus 3:3 6.06 + 0.10 11.69 + 0.24 —0.86 + 0.28 S16 RC (H1) Jones et al. (1999) 

NGC 307 -1.9 8.60 + 0.06 12.02 + 0.24 —0.55 + 0.28 S16 Schw Erwin et al. (2018) 

NGC 821 -4.8 8.22 + 0.19 13.00 + 0.43 -1.27 + 0.46 KH13 GCs PF21 

NGC 1023 -2.6 7.62 + 0.04 12.98 + 0.65 -1.26 + 0.67 KH13 GCs PF21 

NGC 1068 3.0 6.92 + 0.24 12.18 + 0.24 —0.28 + 0.28 S16 RC (Ha) Emsellem et al. (2006) 

NGC 1097 3.3 8.14 + 0.09 12.30 + 0.24 —0.66 x 0.28 vdB16 RC (H1) Ondrechen et al. (1989) 

NGC 1300 4.0 7.88 + 0.30 12.13 + 0.24 —0.68 + 0.28 KH13 RC (H1) Lindblad et al. (1997) 

NGC 1398 2.0 8.03 + 0.08 12.54 + 0.24 —0.53 + 0.28 $16 RC (H1) Moore & Gottesman (1995) 

NGC 1399 -4.6 8.95 + 0.31 12.94 x 0.04 -1.01 x 0.16 S16 GCs Schuberth et al. (2010) 

NGC 1407 -4.5 9.67 x 0.05 13.70 x 0.31 -1.37 x 0.34 S16 GCs PF21 

NGC 2273 0.9 6.93 + 0.04 11.96 x 0.24 —0.60 x 0.28 S16 RC (Hr) Noordermeer et al. (2007) 

NGC 2748 4.0 7.65 + 0.18 11.69 + 0.24 —0.63 x 0.28 KH13 RC (Ho) Erroz-Ferrer et al. (2015) 

NGC2787 -1.0 7.61 + 0.09 12.13 € 0.24 -1.48 + 0.28 S16 RC (H1) Shostak (1987) 

NGC 2960 0.8 7.03 + 0.05 12.43 + 0.24 —0.87 + 0.28 S16 RC (Hr) Sun et al. (2013) 

NGC 2974 -4.3 8.23 + 0.08 12.71 x 0.30 —1.05 x 0.34 S16 GCs PF21 

NGC 3079 6.4 6.40 + 0.05 12.13 x 0.24 —0.82 + 0.28 S16 RC (H1) Sofue et al. (1999) 

NGC3115 -2.9 8.95 + 0.08 13.01 + 0.62 -1.35 + 0.64 KH13 GCs PF21 

NGC 3227 1.5 7.32 + 0.23 12.26 + 0.24 —0.75 + 0.28 S16 RC (H1) Mundell et al. (1995) 

NGC 3245 -2.1 8.38 € 0.11 12.30 + 0.24 —0.97 x 0.28 S16 RC (stars) Zasov et al. (2012) 

NGC 3377 -4.8 8.25 + 0.23 12.22 + 0.34 —0.99 x 0.37 KH13 GCs PF21 

NGC 3607 -3.2 8.14 + 0.15 13.08 + 0.37 —0.96 x 0.40 KH13 GCs PF21 

NGC 3608 -4.8 8.67 + 0.09 13.15 + 0.55 —1.39 + 0.57 KH13 GCs PF21 

NGC 3706 -32 9.77 x 0.06 13.80 + 0.30 -1.81 x 0.34 vdB16 Schw Giiltekin et al. (2014) 

NGC 3783 1.4 7.36 + 0.19 11.76 + 0.24 —0.15 + 0.28 vdB16 RC (H1) García-Barreto et al. (1999) 

NGC 3923 -4.8 9.45 + 0.12 13.44 + 0.15 -1.58 + 0.21 S16 GCs Norris et al. (2012) 

NGC 3998 -2.2 8.93 + 0.05 12.60 + 0.20 —].42 + 0.25 $16 Schw Boardman et al. (2016) 

NGC 4151 1.9 7.81 + 0.08 11.80 + 0.24 —0.33 + 0.28 S16 RC (H1) Mundell et al. (1999) 

NGC 4258 4.0 7.58 x 0.03 12.02 x 0.24 —0.67 x 0.28 S16 RC (Hr) Ponomareva et al. (2016) 

NGC 4303 4.0 6.51 + 0.74 11.76 x 0.24 —0.18 + 0.28 vdB 16 RC (H1) Sofue et al. (1999) 

NGC 4374 -44 8.97 + 0.04 13.69 x 0.57 -1.45 x 0.59 KH13 GCs PF21 

NGC 4388 2.8 6.86 + 0.01 11.96 x 0.24 —0.91 + 0.28 KH13 RC (H^ Woods et al. (1990) 
Veilleux et al. (1999) 

NGC 4459 -1.6 7.84 x 0.08 12.82 + 0.42 —].11 x 0.45 KH13 GCs PF21 

NGC 4472 -4.8 9.40 + 0.04 13.98 + 0.30 —1.69 x 0.34 S16 GCs Cóté et al. (2003) 

NGC 4473 -4.7 7.95 + 0.22 12.88 + 0.51 -1.19 x 0.53 KH13 GCs PF21 

NGC 4486 -4.3 9.81 x 0.06 13.75 x 0.24 —1.40 + 0.28 KH13 GCs PF21 

NGC 4501 3.3 7.30 + 0.08 12.37 + 0.24 —0.60 + 0.28 S16 RC (CO) Nehlig et al. (2016) 

NGC 4526 -1.9 8.65 + 0.12 13.16 + 0.48 —1.17 + 0.50 KH13 GCs PE2I 

NGC 4564 -4.6 7.94 + 0.12 12.88 + 0.78 —1.57 x 0.79 KH13 GCs PF21 

NGC 4594 1.1 8.82 + 0.04 12.73 + 0.24 —0.78 x 0.28 S16 RC (Hr) Bajaja et al. (1984) 

NGC 4649 -4.6 9.67 x 0.10 13.76 x 0.43 -1.43 + 0.46 KH13 GCs PF21 

NGC 4697 -4.5 8.31 +0.11 13.17 + 0.54 -1.29 x 0.56 KH13 GCs PF21 

NGC 4736 2.3 6.83 + 0.12 11.93 x 0.24 —0.76 x 0.28 S16 RC (H1) Speights et al. (2019) 

NGC 4762 -1.8 7.36 x 0.14 12.07 x 0.24 —0.52 +0.28 . Krajnovié et al. (2018) RC (stars) Fisher (1997) 

NGC 4826 2.2 6.19 x 0.13 11.69 x 0.24 —0.28 + 0.28 S16 RC (H1) Braun et al. (1994) 

NGC 4945 6.1 6.13 + 0.18 11.86 x 0.24 —0.70 x 0.28 KH13 RC (H1) Sofue et al. (1999) 

NGC 5328 -4.7 9.67 + 0.16 13.35 + 0.30 —1.31 + 0.34 S16 X-ray Trinchieri et al. (2012) 

NGC 5846 -4.8 9.04 + 0.05 13.85 + 0.45 —1.66 x 0.47 S16 GCs PF21 

NGC 7052 -4.9 8.60 x 0.23 12.91 + 0.30 —0.83 + 0.34 S16 X-ray Memola et al. (2011) 

NGC 7457 -2.7 6.95 + 0.26 12.02 + 0.47 —1.16 + 0.49 KH13 GCs PF21 

NGC 7619 -4.8 9.40 x 0.11 13.82 + 0.37 —1.80 + 0.40 S16 Schw Pu et al. (2010) 


Notes. (1) Galaxy name; (2) morphological T-type from the HyperLEDA database; (3)-(5) black hole mass, stellar mass and star formation efficiency; (6) 
references for BH measurements: KH13-Kormendy & Ho (2013), S16-Saglia et al. (2016), vdB16-van den Bosch (2016); (7) method used for Mp measurements: 
GCs-globular cluster dynamics, RC-rotation curve (using vga, as described in the text), Schw-Schwarzschild model, X-ray- modelling of the X-ray-emitting 
circumgalactic gas; (8) references for Mp measurements. ^ vfat comes from the mean of the two quoted studies. 
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where the mean halo density becomes equal to 200 times the critical 
density of the Universe, that is, My = mgit, 

The first sample that we consider is that of Posti & Fall (2021, 
hereafter PF21), who used measurements of globular cluster radial 
velocities from the SLUGGS Survey (Brodie et al. 2014; Forbes et al. 
2017) to perform a dynamical mass modelling of 25 nearby early-type 
galaxies. The M, and Mp estimates of PF21 rely on the modelling of 
the phase-space distribution of globular clusters using a gravitational 
potential given by the sum of two spherical components, namely a 
Navarro-Frenk-White (NFW; Navarro et al. 1996) dark matter halo 
anda stellar bulge described by a Sersic (1968) profile. The SLUGGS 
data provide astrometric and spectroscopic measurements for several 
tens of clusters per object at distances up to ~ 12 times the galaxy 
effective radii, effectively probing the large-scale galactic dynamics. 
3.6 um images from Spitzer Space Telescope are used to constrain 
the distribution of the stellar components whose mass-to-light ratio is 
a free parameter of the model, along with the halo mass and concen- 
tration. We used the M, and Mp provided by PF21 for a subsample of 
18 of their systems that also have dynamical measurements of Mgy 
collected by Kormendy & Ho (2013) and Saglia et al. (2016). These 
are all based on stellar dynamics, with the exceptions of NGC 4526 
(CO dynamics), NGC 4374 and NGC 4459 (ionised gas dynamics). 

Next, we considered the sample of Terrazas et al. (2017), who have 
combined the data previously collected by Saglia et al. (2016) and 
van den Bosch (2016) to build a dataset of 91 galaxies spanning dif- 
ferent morphological types and with robust estimates for their Mgy 
from stellar, gas or maser dynamics. Stellar masses in Terrazas et al. 
(2017) are derived from extinction-corrected Ks-band photometry 
from 2MASS (Huchra et al. 2012). To determine the Mp of these 
objects, we have searched the literature for spatially resolved obser- 
vations of large-scale kinematic tracers, preferably H 1 emission-line 
data but also Ha or CO data. In cases where we judged the velocity 
field of such tracer to be regular enough so that a rotation curve could 
be determined reliably, we used vga, the rotational speed in the flat 
part of the rotation curve, as a proxy for Mp. The details of the vfat- 
to- My calibration are presented in Appendix A. When gas kinematics 
was not available, which is typically the case in earlier galaxy types, 
we used dynamical estimates based preferably on globular cluster 
radial velocities or, in the absence of these measurements, on mod- 
elling of the X-ray emitting halo gas, or on Schwarzschild models for 
the stellar component. Following this approach, we included a total 
of 30 galaxies from the Terrazas et al. (2017) sample. Finally, we 
complemented our dataset with three additional objects from the re- 
cent sample of de Nicola et al. (2019), and four spiral galaxies from 
Kormendy & Ho (2013), following the same procedure described 
above to estimate the halo masses. 

The main properties of the resulting sample of 55 galaxies are 
listed in Table 1, along with the references related to the Mgy and 
My measurements. The vast majority of galaxies in our sample have 
Mp determined with the rotation curve method (27, of which 21 
from Hı data, three from Ha, one from CO and two from stellar 
kinematics) or from globular cluster dynamics (22). Schwarzschild 
models and the X-ray method are used only for four and two objects, 
respectively. We stress that the Mp measurements used in this work 
have not been derived from observations alone, since none of the data 
extend to the halo virial radii. Instead, they are based on theoretical 
constraints on the properties of DM halos in ACDM cosmological 
models, and in particular on the correlation between halo mass and 
concentration of Dutton & Macció (2014) which our rotation curve 
method and the measurements of PF21 rely on. Aside from the 18 
early types from PF21, where M, is determined dynamically, we have 
homogenised the stellar mass measurements in our sample using a 
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common M, / Ly ratio of 0.6 (McGaugh & Schombert 2014). Other 
choices for Mx / Ly are possible, but have little impact on our results. 

As most galaxies have Mp determined either from rotation curves 
or from globular cluster dynamics, it would be important to quantify 
how these two methods compare to each other. This is not trivial, 
given that globular clusters are more abundant in early type galax- 
ies where cold gas is typically scarce. In our sample, NGC 2974 
is the only galaxy for which both measurements are available: us- 
ing our eq. (A1) with vfat =355 + 60 km s^! (from the H 1 study of 
Kim et al. 1988) gives log}g(My/ Mo) = 12.68 + 0.27, which per- 
fectly agrees with the value found by PF21 and reported in Table 1, 
logj9(Mp/Mo) = 12.71 x 0.30. This result, even if it is obtained 
for a single object, is reassuring and preparatory for the rest of the 
analysis. 


3 OBSERVED SCALING RELATIONS 


We now focus on the relations between fx, Mpy, and My for the 
systems in our sample. In Fig. 1 we show the distribution of our 
galaxies in the (Mpg. Mp) space (left panel), in the (f£, Mp) space 
(central panel) and in the (f;,Mpy) space (right panel). It is evident 
that, in the mass range considered, each quantity appears to be related 
to the others. On average, the star formation efficiency decreases with 
increasing halo or BH masses, while Mp and Mpy are positively 
correlated with each other. 

The correlation between Mpy and Mp does not come as a surprise, 
as it was firstly discovered by Ferrarese (2002). While we discuss 
this relation further in Section 5.1, we highlight here a couple of 
interesting features. First, the trend shown by the data at high Mp 
seems to break down around Mp of a few x10? Mo (or Mpy of 
107 — 108 Mo), steepening below this threshold mass. This feature 
may be simply produced by an increase in the scatter of the relation 
in the low Mp regime, possibly coupled with low-number statistics, 
but we will later show that a break in the Mpy-M, relation arises 
‘naturally’ in our galaxy evolution model as resulting from a change 
in the mode by which BHs accrete their gas. Second, the observed 
correlation is surprisingly tight. In Section 3.2 we show that, in our 
sample, the strength of the correlation between Mpy and Mp and 
its intrinsic scatter are comparable to those of a very well studied 
scaling law, the Mpy-o relation. 

The relation between f, and Mhalo, shown in the central panel 
of Fig. 1, is also not surprising. Galaxies are well known to fol- 
low a stellar-to-halo mass relation (SHMR, e.g. Moster et al. 2013; 
Behroozi et al. 2013) according to which f, peaks at the value of 
~ 0.2 for halo masses of ~ 10!* Mo and decreases rapidly at lower 
and higher Mp, possibly because star formation is made inefficient 
by stellar and AGN feedback, respectively. As the minimum Mhp in 
our data is close to this peak value, we only sample the high-mass, 
descending portion of the SHMR. Interestingly, some low-Mp galax- 
les in our sample have f, compatible with 1, meaning that they have 
been able to convert all their (theoretically) available baryons into 
stars. The same result was found by Posti et al. (20192) for a sample 
of high-mass spirals from the SPARC dataset (Lelli et al. 2016). 

The rightmost panel of Fig. 1 shows the anti-correlation between 
fx and Mpy. As noticed for the Mgy-My relation, also here the 
trend seems to change around Mpy ~ 107? Mo: galaxies hosting 
lower-mass BHs have approximately constant fx, while at larger 
Mpu the star formation efficiency gets progressively reduced. As 
discussed before, the low-mass trend could also be due to an increase 
in the scatter, although here the break appears more evident. As this 
relation is less tight than that between fx and Mp, one may conclude 
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Figure 1. Mgy-My relation (leftmost panel) fx-Mp relation (central panel), and fx-Mpgy relation (rightmost panel) in our galaxy sample. Markers in each 
panel are colour-coded according to the third dimension in the (My, Mpn. fx) space. Circles (squares) are used for systems with Hubble morphological T-type 
< 0 (> 0), corresponding to early (late) galaxy types. Star markers are used for galaxies featuring pseudo-bulges in Kormendy & Ho (2013). As a reference, the 


locations of the Milky Way (MW) and M31 are shown. 
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Figure 2. Top panel: residuals in the f, - My relation vs residuals in the Mpy- 
Mh relation for our galaxy sample. At fixed Mj, there is weak tendency for 
galaxies hosting heavier BHs to show larger f, (Pearson coefficient of 0.1). 
The colour palette shows the Mp of each galaxy. Bottom panel: residuals in 
the fx-Mpgu relation vs those in the Mp-Mpgy relation. At fixed Mpg, there is 
a strong evidence for galaxies living in lighter halos to have larger f (Pearson 
coefficient of —0.6). The colour palette shows the Mgy of each galaxy. The 
markers used in both panels are the same as in Fig. 1. 


that mass, and not BH feedback, is the main quenching driver in high- 
mass galaxies (e.g. Bundy et al. 2008; Peng et al. 2010b; Geha et al. 
2012; Dubois et al. 2013). The results from our modelling (Section 
4), however, suggest that this is not the case. 


In Fig. 1 we have used different markers for galaxies of different 
morphology: circles show earlier galaxy types (T < 0), squares show 
later galaxy types (T > 0), while star markers highlight the spirals 
hosting pseudo-bulges as listed by Kormendy & Ho (2013). While 
early-type systems populate preferentially the high Mp regime and 
late-types appear only at lower Mj, the galaxy population as a whole 
seems to distribute along a well defined sequence in the Mgy-Mp- f 
space. Discs hosting pseudo-bulges are no exception, as they do not 
occupy a preferential position in any of the three relations presented. 
We stress that the overlap between different galaxy types in the fx- 
Mp, plane (central panel of Fig. 1) is not in tension with the results 
of PF21, who found distinct SHMR for early and late type systems: 
the split between morphological types becomes evident only in the 
f«-M, space (not shown here), due to the fact that the most massive 
early and late-type galaxies, albeit having halo masses that differ by 
more than an order of magnitude, show only a factor 2—3 difference 
in M, (see Fig.3 in PF21). 


The color-coding used in Fig. 1 is useful to highlight another in- 
teresting feature of our data. The first and the third panel of Fig. 1 
suggest that, at fixed Mpy, systems that have been more (less) effi- 
cient at forming stars are those that live in lower (higher) mass halos. 
Conversely, no significant trend can be seen at fixed Mp (first and 
second panel). The features discussed can be better appreciated in 
Fig. 2. Here, we have determined empirically the mean trends of fx 
and Mgg as a function of Mp, and plotted the residuals with respect 
to these trends against each other in the top panel of Fig.2. The 
clear lack of an anti-correlation between these residuals is surpris- 
ing, as one might expect that galaxies hosting larger BHs are also 
those where AGN feedback, and therefore star formation quenching, 
is more effective. In contrast, there is instead weak evidence for a 
positive correlation (Pearson coefficient of 0.1). Similarly, the resid- 
uals with respect to the mean f,-Mpy and Mp-Mpgqų relations are 
compared in the bottom panel of Fig. 2. The anti-correlation between 
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Figure 3. Two different 3D views for the distribution of our galaxy sample 
in the (My, Mpg. fx) space. The markers used in both panels are the same as 
in Fig. 1. The best-fit plane from eq. (1) is shown in dark grey. The 3D curve 
passing through the data shows the fiducial theoretical model from Section 4. 
To better enhance the 3D perspective, galaxies are colour-coded by their fx 
and 2D projections on the (Mi, Mpgg) plane are shown. 


fx and Mj at fixed BH mass becomes now very evident (Pearson co- 
efficient of —0.6). In Section 4.2 we show that uncorrelated scatter in 
the parameters of our evolution model can explain the trends shown 
in Fig. 2 remarkably well. In particular, system-to-system variability 
in BH feedback efficiency is a viable explanation for the observed 
anti-correlation between f, and Mj at fixed Mpy. 


3.1 The 3D Mgy-M,-f, relation 


We now analyse the distribution of our galaxy sample in the 3D 
(fx, Mgp. My) space. The eingenvalues of the covariance matrix as- 
sociated to our data are approximately in a 26:3:1 ratio. The presence 
of a substantially larger eigenvalue indicates that there is a clear ten- 
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dency for the data to distribute along a unique direction in the 3D 
space, given by the associated eigenvector. This is a simple geomet- 
rical confirmation of the fact that Mgy, Mp and fx are all related to 
each other, as previously discussed. It also suggests that, in the mass 
range studied here, measuring either Mpy, Mx or Mp in a galaxy 
suffices to determine the other two quantities, with good approxima- 
tion. 

The 3:1 ratio of the two lower eigenvalues indicate that a 
(marginally) more refined representation for our data is given by 
a plane in the 3D space considered. We used the LrsFrr Python pack- 
age from Cappellari et al. (2013) to fit our data with the following 
parametric form: 


M, M 
log( fx) = a logo (a) + b logio (=x) +c. 


(1) 
The code uses a least-squares fitting algorithm which allows for 
intrinsic scatter and errors in all coordinates. The best-fit solution 
is found for a = —0.77 + 0.15, b = 0.13 + 0.09, c = 7.70 + 0.05, 
and is consistent with no intrinsic scatter. The parameters found 
are indicative of the same trends discussed above, i.e., f, depends 
strongly on My at a fixed Mgy, but very weakly (or not at all) on 
Mpu at a fixed Mp. The best-fit parameters of eq. (1) remain within 
the quoted uncertainties if we fit only the data with Mp > 10? Mo 
(that is, we exclude the data points below the break). 

Figure 3 offers two representative 3D views of the data, along 
with the best-fit plane and the theoretical model that we build in 
Section 4. The 3D views clearly highlight how the data distribute 
preferentially along a one-dimensional sequence, although a plane 
can capture the distribution of their scatter. We stress, though, that 
the plane described by eq. (1) is valid exclusively in the range of 
masses spanned by the data, and in particular it will not hold at 
My < 5x 10!! Mo, a regime where f, is known to decrease, rather 
than increase as eq. (1) would suggest. 


32 Relating Mgg to the properties of the stellar component 


The importance of characterising the relations between Mpy and 
other observable galaxy properties is twofold. On the one hand, 
these relations provide fundamental clues to constrain the physics 
of BH growth and to clarify the role of AGN feedback in galaxy 
evolution. On the other hand, the existence of tight scaling relations 
offer convenient ways to determine the BH masses using more easily 
observable quantities as proxies. While the primary focus of this 
study is the relation between Mpg and global galaxy properties such 
as My and fx, in this Section we briefly discuss how BH masses 
relate to some of the properties of the stellar component in our 
galaxy sample. 

Fig. 4 shows the relations between Mgg and four different stellar 
properties: the total stellar mass Mx, the stellar disc mass Mgisc, 
the stellar bulge mass Mpulge, and the mean stellar velocity disper- 
sion within the effective radius ce. The Mpg — Mp relation, already 
discussed in Section 3, is also shown as a reference. The og measure- 
ments adopted here are from de Nicola et al. (2019) and Kormendy 
& Ho (2013), while bulge fractions are mostly based on the 2D- 
decompositions of 3.6 um images from the the Spitzer Survey of 
Stellar Structure in Galaxies (Sheth et al. 2010) via the GALFiT pack- 
age (Peng et al. 2002, 20102) or, when these were not available, 
on the kinematic decomposition reported in Fall & Romanowsky 
(2018). Each panel in Fig. 4 reports the Spearman (S) and Pearson 
(P) correlation coefficients of the quantity pair analysed, and the 
best-fit parameters of linear (in log-space) regression determined 
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Figure 4. Relations between Mpy and other galaxy properties in our sample: halo mass Mp (top-left panel), stellar mass M. (top-middle panel), stellar disc 
mass Maise (top-right panel), stellar bulge mass Mpulge (bottom-left panel) and mean stellar velocity dispersion within the effective radius o. (bottom-right 
panel). Symbols are the same as in Fig. 1, colour-coded by the Hubble morphological type of each galaxy. In the top-left corner of all panels we list the Spearman 
(S) and Pearson (P) correlation coefficients of the data, followed by the intercept (a), slope (b), vertical and perpendicular intrinsic scatter (Avert and Aperp) 
of the best-fit linear regression determined with LrsFrr (Cappellari et al. 2013). The best-fit relation is shown as a light-blue solid line, with the filled region 
corresponding to +Avert around it. Dotted and dashed lines in the various panels show some of the scaling relations previously reported in the literature and are 


discussed in Section 5.1. 


with LrsFrr, including the intrinsic scatter in the vertical and per- 
pendicular directions, Avert and Aperp. Galaxies in Fig. 4 are colour 
coded by their Hubble morphological T-type, taken from the Hyper- 
LEDA database. Fig. 4 also shows a representative selection of the 
best-fit relations previously determined in the literature, which will 
be discussed in more detail in Section 5.1. 


The correlation coefficients that we find confirm the widely ac- 
cepted scenario according to which i) the stellar property that best 
correlates with Mpy is Ge; ii) bulges correlates with Mpy much bet- 
ter than discs. However, in our sample, the strength of the correlation 
between Mpy and Mj (correlation coefficients of ~ 0.83) is com- 
parable to that between Mpy and oe (~ 0.86). We also find a lower 
Avert in the Mpy—M,y relation than in the Mpg — ce relation, although 
this is probably due to the large uncertainties associated with our Mp 
estimates. All considered, our results indicate that both ce and Mp 
seems to provide a similar accuracy when used as proxies for the BH 
mass, which is quite remarkable given the completely different scales 
involved. We stress, though, that the perpendicular intrinsic scatter 
Aperp of the Mgy-ce is about half that of the Mgy-Mp, indicating 
that the former is more ‘fundamental’ than the latter. The relation 
between c and Mp (not presented here) shows a trend similar to that 
between Mpy and Mp, with a hint of a break visible around c of 
130-160 km s^! . As discussed by PF21, a correlation between these 
two quantities is expected as resulting from the combination of the 
SHMR, which links Mp to Mx, and the Faber & Jackson (1976) and 
Tully & Fisher (1977) relations, which relate Mstar to a characteristic 
velocity of the stellar component (that is, Ce for early galaxy types). 


Another interesting feature shown by the data is that the Mpy- 
Mbulge relation shows similar strength and intrinsic scatter as the 
Mgg — Mx relation. This means that, in our sample, the bulge mass 
and the total stellar mass are equally-good proxies for Mgy, which is 
a puzzling result considering the poor correlation between Mpy and 
Maisc. The top-right panel of Fig. 4 may give a clue to the solution 
of this puzzle, indicating that earlier and later galaxy types follow 
somewhat parallel sequences in the Mgy-Maisc plane. This results in 
an overall poor correlation between BH and disc masses, but when 
bulges are included in the total M, budget the sequence of early 
types approaches the other one, increasing the overall correlation 
strength. A similar result was also found by Davis et al. (2018), who 
stressed the importance of distinguishing between early- and late- 
galaxy types in the study of the Mgy — M, (see also Section 5.1). We 
note, however, that relative to a unbiased, volume-limited sample, 
massive discs are over-represented (shallower mass function) in our 
sample while massive spheroids are under-represented (steeper mass 
function). This limitation, together with the fact that bulge/disc frac- 
tion measurements are available only for 70% (39/55) of it, restrain 
us from investigating these features further in this work. 


Finally, we stress that the correlation coefficients and best-fit pa- 
rameters determined in this Section do not vary significantly when 
galaxies hosting pseudo-bulges are removed from the analysis. The 
main difference is visible in the Mpy — Mp relation, whose slope 
decreases when these low-mass systems are removed. However, this 
would be readily explained if the intrinsic shape of the Mgy — Mp 
relation was not a simple power-law, but its slope increases at lower 


MNRAS 000, 1—20 (2021) 


8 A. Marasco et al. 


masses as we suggest in the Section below. Similarly, moderate 
(^ 30%) variations in the mass-to-light ratios of discs and bulges 
have little impact on the results presented in this Section. 


4 A SIMPLE GALAXY EVOLUTION MODEL 


In this Section we investigate the physical origin of the trends pre- 
sented in Section 3 using a simple equilibrium model for galaxy 
evolution in the ACDM framework!. Our model is largely inspired 
by the work of B17 and is based on a commonly accepted framework 
where galaxy halos smoothly accrete dark matter and gas at a cos- 
mological rate, having their stellar and black hole build-up regulated 
both by the cooling of the available gas reservoir and by stellar/AGN 
feedback. We provide a full description of our model in Appendix B, 
while below we give a brief overview - sufficiently detailed to follow 
the rest of this study - of its main ingredients. 

The model follows the evolution of BH, stellar and dark matter 
masses of a galaxy given some initial conditions, namely the initial 
‘seed’ masses Mpy seed and Mp seed, and the seeding cosmologi- 
cal time fseeq. At each time-step, the mass growths of the various 
components are determined via the sequence of physical processes 
illustrated in Fig. 5. These include standard recipes for cosmological 
accretion of dark matter and gas onto halos, gas radiative cooling, 
star formation, BH accretion, and simplified treatments for stellar 
and AGN feedback. The galactic gas reservoir, intended as the sum 
of interstellar and circumgalactic media (ISM and CGM), is de- 
scribed as a single component following an equation of state with 
y = 4/3 (as in B17), with a temperature equal to the halo virial 
temperature, a mass given by? fo My — My — Mpun, and a pristine 
metallicity (chemical enrichment is not treated in our model, but a 
case with a higher metallicity is discussed in Section 4.2). At each 
time-step, the gas cooling rate is balanced by the rates at which stars 
form, the BH grows and gas is returned to the initial reservoir be- 
cause of stellar mass losses and feedback (‘equilibrium’ model). In 
our implementation, feedback from star formation has two effects: it 
drives galaxy-scale outflows which instantly return gas to the initial 
reservoir, with a mass-loading 6 equal to (My/ Mei) ^ (Merit and 
a being free parameters), and it acts as a regulator of the gas density 
close to the BH, significantly reducing its accretion rate for large 8 
(e.g. Hopkins et al. 2021). AGN feedback is treated as a continuous 
accumulation of energy deposited by the BH onto the gas reservoir 
at a rate cc epMpg, er being the BH feedback efficiency, another 
free parameter. The total energy released by the BH, Egy, must not 
exceed the gravitational binding energy of the cooling gas, Eggo: at 
time-steps when this occurs, the BH accretion rate is reduced so that 
the condition above is satisfied, and star formation is manually turned 
off. This form of 'preventative' AGN feedback is the only quench- 
ing channel that we provide in our galaxy evolution framework. The 
main parameters regulating our model are summarised in Table 2, 
while the caption of Fig. 5 points to the relevant equations defined in 
Appendix B. 

We stress that the final goal of this modelling exercise is to offer 
a simple but cosmologically motivated interpretation for the ob- 
served scaling relations between BH, stellar and halo masses pre- 
sented in Fig. 1 and 3. While we acknowledge that Mpy correlates 
more strongly with Mpulge than with Maise (Section 3.2), here we 


! We assume a Planck Collaboration et al. (2013) cosmology (Qm,0 = 0.315, 
Ho = 67.3 kms"! Mpc™!), in particular f; = Qp/Qe = 0.188. 

? Hence the total mass of baryons within the halo is always equal to fj, My, 
which implies that the gas accreted onto the halo never leaves the system 
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do not offer separate treatments for the growth via smooth accre- 
tion or via episodic mergers and treat the galaxy stellar component 
as a whole. This simplifies our approach significantly, bypassing 
the need for prescriptions regulating the formation of bulges which 
would introduce additional complexity to the model. Clearly, the use 
of such simplification implies that we will not get physical insights 
on the co-evolution between bulges and BHs. For instance, we cannot 
distinguish a scenario where the Mgy-Mpulge relation results from 
the physics of stellar and AGN feedback, which separately produce 
Mpy « c^ and Mbulge * ot scalings (Power et al. 2011; King & 
Nealon 2021), from a scenario where the Mgy-Mpbuige arises statisti- 
cally from the hierarchical merging of galaxies with initially uncor- 
related Mgy and Mpulge (Peng 2007; Jahnke & Macció 2011). We 
leave the answer to this topic to more dedicated theoretical studies, 
though our simple model still provides useful insights once com- 
pared to the presented correlations between Mpy and global galaxy 
properties. 


4.1 A fiducial model 


Table 2 shows a summary of the main parameters of our equilibrium 
model. By experimenting with different parameter values, we have 
found a ‘fiducial’ set which produces a model in excellent agreement 
with our data at z 2 0. The fiducial parameter set is reported in the 
bottom portion of Table 2 and the comparison with the data is pre- 
sented in Figure 6. Clearly, our model predicts a Mgy-Mp- fx relation 
(dashed lines in Fig. 6) which passes right through the data. In par- 
ticular, the model predicts a break in the relations at Mpy = 107 Mo 
which is visible in our data as well, as discussed in Section 3. In the 
central panel of Fig. 6 we also show the SHMR determined by Moster 
et al. (2013) from abundance matching prediction (green-shaded re- 
gion), and the dynamical measurements of Posti et al. (20192) for 
the spirals from the SPARC dataset (Lelli et al. 2016), which allows 
us to extend the dynamical range of the fx-Mp plot down to lower 
masses. However, measurements for Mpy in the SPARC sample are 
not available. Our model performs very well along the entire mass 
range. Interestingly, the measurements from SPARC indicate that the 
scatter in the SHMR increases at lower Mp, a point which we will 
return to in Section 4.2. A 3D view of the fiducial model was already 
offered in Fig. 3. 

While our model seems to correctly predict the relations between 
stellar, halo and BH masses, which are quantities integrated over cos- 
mic time, one may question its performance when it comes to more 
‘instantaneous’ properties, like the star formation rates of present-day 
galaxies. In our model, the SFR of a system is abruptly shut down 
when the BH energy output becomes equal to the binding energy of 
gas within the halo cooling radius. We show below that this happens 
earlier in more massive halos, while in less massive galaxies this may 
only occur at £ > tyypple- This produces a segregation between low- 
mass systems that are normally star forming at z=0, and high-mass 
galaxies that are quenched. Figure 7 shows the relation between the 
SFR, averaged over the latest 100 Myr of evolution, and the stellar 
mass of z = 0 systems in our fiducial model, and compares it with 
the approximate observed distribution from SDSS measurements by 
Renzini & Peng (2015). Our model correctly reproduces the slope 
(^ 1) of the main sequence of star formation (e.g. Noeske et al. 
2007; Popesso et al. 2019), albeit with a slightly lower normalisation 
(by ~ 0.2 dex), as well as the approximate M, at which quenched 
galaxies begin to dominate the SFR-M, distribution (that is, around 
10105 Mo). Clearly, the transition between these two regimes is 
supposed to be gradual, whereas in our simplified treatment of BH 
feedback it is abrupt. We discuss this further in Section 5. Going 
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Figure 5. Sketch showing the sequence of physical processes regulating the mass build-up in our galaxy evolution framework within a time interval [t, t + At]. 
(i) At a time ¢ the galaxy is fully described by its My, Mpu and Mx. Its gas content Mgas (ISM+CGM) is such that the total baryonic mass is fj My. The 
gas temperature is set to the halo virial temperature. (ii) Dark matter and gas accrete onto the halo at a cosmological rate given by eq. (B1). (iii) A fraction of 
Mgas cools, instantaneously providing fuel for star formation and BH accretion. The cooling mass depends on the variation of the gas mass enclosed within the 
cooling radius (7991, eq. B2) in the time interval considered. (iv) A first estimate for the mass accreting onto the BH, AMgpu, is computed via eq. (B7) and (B8). 
A first estimate for the mass of the newly formed stars AM, is computed via the equilibrium model, eq. (B3). The energy injected by the AGN feedback into the 
gas reservoir is increased by an amount AEgg given by eq.(B9). (v): Egu is compared to the gas gravitational binding energy within the cooling radius, Ecool 
(eq. B10). If Egy < Eecool, the BH mass is increased by AMgg and the stellar mass by AM, , otherwise AMgy is reduced so that Egy = Ecoo) and AM, is set 


to zero. 


Table 2. Main parameters of our theoretical model of galaxy evolution. The values listed for the free parameters are those of our fiducial model discussed in 


Section 4.1. 

Parameter description value status 
fsesdi min earliest seed injection time 0.3 Gyr (z 213.9) fixed 
teed; max latest seed injection time 4.0Gyr(z 21.62) fixed 
Mh seed seed halo mass 1019 Mo fixed 
M, seed seed stellar mass 10? x MBH, seed fixed 

Z gas gas metallicity pristine fixed 
R recycled gas fraction due to stellar mass losses 0.3 fixed 
MBH. seed seed black hole mass 5x 10*Mo free 
DBH,0 normalisation of the gas density near the BH, eq. (B8) 0.35 cm? free 
€f BH feedback efficiency, eq. (B9) 1.0 x 107? free 
Mcerit,o present-day critical halo mass, eq. (B5) 2.5 x 10!! Mo free 
Qa slope of the mass-loading due to stellar feedback, eq. (B4) 1.7 free 


back to Fig. 6, we have used blue and red dashed-lines to show the 
star-forming and quenched galaxy populations, respectively. In our 
fiducial model, at z = O the transition occurs at My, = 10? Mo, so 
that most of the systems studied in this work are supposed to be 
quenched. While this may not be the case for individual galaxies, our 
results hold for the galaxy population as a whole. 


Further insights on the evolutionary scenarios predicted by our 
fiducial model are provided by Figure 8, where we show the dark 
matter, stellar and BH mass build-up as a function of z for systems 
of different present-day Mp. By construction, the halo growth (panel 
a in Fig. 8) proceeds undisturbed at all z and for all systems. The 
BH growth (panel c), instead, is more interesting. At high z the BH 
accretion proceeds slowly, both because seeds are still small and 
because supernova-driven outflows are very efficient at lowering the 
central gas density. As stellar feedback becomes less efficient, this 
early phase is then followed by a very rapid growth, which ends 
abruptly only when Epy = Egoo}. Then, the accretion proceeds at a 


more gradual pace in a self-regulating mode driven by the continuous 
balance between Egy and E.o,). This balance defines univocally the 
slope of the Mgg-Mj relation for Mp = 10? Mo, whose excellent 
agreement with the data is one of the main success of our simple 
theoretical model. The z at which the transition to the self-regulating 
accretion mode occurs is a function of the galaxy mass: high-mass 
systems enter the self-regulating phase earlier, whereas galaxies with 
My(z 20) € 10!2Mo (corresponding to Mpy(z=0) € 107 Mo) 
are still in the rapidly accreting mode. This agrees well with the 
picture of anti-hierarchical growth of BHs emerging from the study 
of AGN luminosity functions (e.g. Fig. 8 in Marconi et al. 2004), and 
produces the break in the scaling relations visible around BH masses 
of 107-10? Mo (see Fig. 1). Also, this transition between the rapidly 
accreting phase and the self-regulating phase signals the end of the 
stellar mass build up (panel b). As the transition occurs earlier in more 
massive halos, the redshift at which quenching occurs increases for 
increasing M, (z 20), which is again in agreement with the picture 
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Figure 6. Fiducial model at z = 0 (thick dashed curves) vs. observational data from this work. Markers are the same as in Fig. 1, but colour-coded by the Hubble 
morphological type of each galaxy. Left panel: Mgy vs My. Central panel: fx vs My, triangles are measurements from Posti et al. (20192) using SPARC spirals 
(Lelli et al. 2016), the green-shaded region shows the abundance matching prediction from Moster et al. (2013). Right panel: fx vs Mpu. In all panels, the blue- 
(red-) dashed curve show model galaxies that are star-forming (quenched) at z = 0. 
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Figure 7. SFR averaged over the last 100 Myr vs M, at z = 0. The blue 
dashed line and the red triangles show galaxies from our fiducial model that 
are star-forming or quenched, respectively. The SFR of quenched system 
is set to an arbitrarily low value. The blue (red) shaded region shows the 
approximate distribution for star-forming (quenched) systems in SDSS from 
Renzini & Peng (2015). 


of anti-hierarchical growth of galaxies (or ‘downsizing’, e.g. Cowie 
et al. 1996; Eyles et al. 2005; Fontanot et al. 2009). However, in 
our model the BHs of quenched galaxies still accrete gas because 
My and E95; keep growing with time, requiring increasingly higher 
Epu for the balancing. Thus the BH-to-stellar mass ratio increases 
with time in most massive halos (panel d), which seems to be at odds 
with predictions from other theoretical models (e.g. Lamastra et al. 
2010) and with observational data (e.g. Merloni et al. 2010), although 
observational bias at high-z may play an important role (Lauer et al. 
2007). We discuss further the limitations of our modelling approach 
in Section 5.3. 

The reader may have noticed that we have avoided any statistical 
approach to optimise the model parameters to the observed data. In 
fact, the fiducial values reported in Table 2 are only indicative, and 
we do not exclude that a better match could be obtained by tuning 
the parameters further or even using very different parameter values. 
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The spirit of our approach is to offer a proof of concept that a simple 
equilibrium model in ACDM framework is well suited for describing 
many of the observed trends, rather than offering precise estimates 
for some parameter values. 


4.2 Interpreting the observed scatter 


One of the most interesting features of our data lies in their scatter, 
and in particular in the correlations between the residuals shown in 
Fig. 2. Inthis Section we explore the possibility that such correlations 
arise ‘naturally’ in our model because of random fluctuations in the 
model parameters. 

As a first step, it is instructive to assess how the Mpu-Mp- fx re- 
lation predicted at z = 0 responds to variations in the single model 
parameters. This is shown in Figure 9, where we compare our fidu- 
cial model with a series of other realisations obtained by varying 
one by one the various parameters while keeping the others fixed to 
their fiducial values. This exercise reveals a number of interesting 
features. First, the high-mass slope of the Mgg-Miy relation remains 
the same in all our experiments, indicating that it is a strong predic- 
tion of the self-regulated accretion mechanism described above. Its 
normalisation instead depends entirely on er (panel c in Fig. 9): as 
expected, higher BH masses are reached for lower feedback efficien- 
cies, and vice versa. Surprisingly, er has no impact on the SHMR. 
As we discuss in Section 5.3, this is due to the primary role of BH 
in processing the gas accreted from the halo and can be seen as a 
consequence of our oversimplified treatment of star formation pro- 
cesses. In practice, ef regulates the fraction of Moo that feeds the 
BH in the self-regulating phase, but (by construction) the fraction 
that does not accrete onto the BH does not form stars either, leav- 
ing no impact on the SHMR. The role of MBH,seea and ngy,o is 
to set the BH accretion rate via eq. (B6) during the phase that pre- 
cedes the self-regulating mode, effectively shortening or extending 
the self-regulating sequence in the Mgy — My plane and affecting the 
post-peak part of the SHMR (panels a and b). We note that Mp seed 
and ngy,0 are completely degenerate in the mass range considered, 
which is another reason for considering the values quoted in Table 
2 as indicative only. Finally, the parameters regulating the stellar 
feedback efficiency, a and Merit,o, play a major role in setting the 
slope of the SHMR at all masses and its normalisation at low masses, 
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Figure 8. Dark matter (a), stellar (b) and black hole (c) mass build-up as a function of redshift in our fiducial model, colour-coded by the halo mass reached at 


z = 0. Panel (d) shows the evolution of the BH-to-stellar mass ratio. 


respectively (panels e and d). Thus, in our models, stellar feedback 
participates in regulating the star formation efficiency at all masses 
(via æ), and not only in the low-Mj regime. 

The panel-set f in Fig.9 shows the effect of altering some of the 
main assumptions of the model, and deserves a more in-depth discus- 
sion. The coefficient of the equation of state for the gas accreting onto 
the BH has an influence on the slope of the SHMR in the high-mass 
regime: using y = 1 in eq. (B7) (isothermal gas) leads to a steeper 
slope (purple-dashed curves) that is not compatible with the data 
and that cannot be easily changed by varying the other parameters 
of the model. Removing completely the BH feedback from our fidu- 
cial model (that is, assuming er = 0) leads to overmassive BHs in 
all galaxies, as shown by the red-dashed curves in the panels. Al- 
lowing the BH growth to self-regulate without suppressing also the 
star formation in systems with Egy > Eeool leads to a much flatter 
high-mass slope for the SHMR (blue-dashed curves), which means 
that star formation quenching driven by BH feedback is an important 
ingredient in our model. Finally, the green-dashed curves show the 
effect of increasing the gas metallicity from pristine to 0.1 Zo. This 
promotes gas cooling and effectively enhances the star formation ef- 
ficiency, increasing the maximum f, achievable without changing 
the peak mass. This would give a better agreement the high- fx values 
measured for the most massive spirals of the SPARC sample (Posti 
etal. 20192), and for some of the systems in our sample as well. While 
modelling the metallicity evolution of the galactic gas reservoir goes 
beyond the purpose of the present study, we stress that the peak fx 
reached by our fiducial model may well be underestimated given the 
pristine composition assumed for the gas accreted at all redshift. 

The above experiments offer an interpretation of the scatter seen 
in the data. From panel-set c in Fig. 9, it is clear that a scatter in 
er produces an anti-correlation between f, and Mp at a fixed Mpy. 
This happens because, as eg scatters from lower to higher values, 
galaxies hosting BH of similar mass occupy more massive halos 
(that is, generated from earlier seeds), which necessarily correspond 
to lower values of f, as the SHMR is not sensitive to er. Fluctuations 
in Miro, instead, offer a way to obtain both a positive correlation 
between fx and Mpu ata fixed Mp (but only in the low-Mp regime), 
and a higher scatter in the low-M,, part of the SHMR. 

These considerations can be better visualised in Fig. 10, which 
shows a stochastic realisation of our fiducial model obtained by in- 
troducing Gaussian fluctuations, with standard deviation of 0.4, in 
log(er), log(Merit,9/ Mo) and a. The color-coding follows that of 
Fig. 1 and helps one recognise, within the model, the same features 
seen in the data and discussed in Section 3. By fitting eq. (1) to this 
synthetic dataset, after excluding systems with Mp <5 x 10!! Mo or 
Mpu < 106 Mo (grey circles in Fig. 1) which do not have a counter- 


part in the observed sample, we find a 2 —0.68, b 20.05 and cz 7.21. 
These value are in excellent agreement with those found for our 
Observed dataset in Section 3.1. In Fig. 11 we show the relations be- 
tween the residuals calculated at fixed Mp and Mpy in our scattered 
model, computed following the same procedure used for the data in 
Fig. 2. The comparison between Fig. 2 and Fig. 10 shows remarkable 
similarities between the data and the model (similar Pearson coeffi- 
cients). Including also fluctuations in MgH,seea and ppp, o does not 
alter the trends seen in the residual, although the agreement with 
the data worsens. Although qualitative, these considerations indicate 
that a simple, uncorrelated scatter in the model parameters may well 
explain the features observed the data residuals. 

Once again, the experiments presented in this section must be 
considered as a proof of concept that the data can be well reproduced 
by models like the one presented in this work. The physical origin 
for the scatter in the model parameters is a matter of great relevance, 
but its investigation goes beyond the purpose of this study. 


5 DISCUSSION 
5.1 Comparison with other works 


The positive, tight correlation between Mp and Mpy was first no- 
ticed by Ferrarese (2002) as a consequence of the correlation between 
rotational velocity and central velocity dispersion in galaxies (origi- 
nally discovered by Whitmore et al. 1979), the former being a proxy 
for Mp and the latter for Mgy. Later on, Bandara et al. (2009) pro- 
vided a more direct evidence for a correlation between Mpy and the 
total dynamical mass of the galaxy Mayn by focusing on a sample 
of early type systems for which Mayn were determined via gravita- 
tional lensing models. A similar Mgy-Mayn relation was derived by 
Krumpe et al. (2015) by studying the X-ray luminosity dependence 
on the clustering strength of low-z AGNs. Interestingly, the original 
Mgg- Mj relation determined by Ferrarese (2002), shown as a black 
dotted line in the top-left panel of Fig. 4, is in a remarkably good 
agreement with our data, in spite of the different methods used to 
determine BH masses and the different v fat — Mp calibration adopted, 
and with the results of Bandara et al. (2009) and Krumpe et al. (2015). 

The findings of Ferrarese (2002) were very surprising as they 
indicated that dark matter alone could engineer the BH growth with- 
out passing through the complex baryonic physics associated with 
BH accretion. Kormendy & Bender (2011) and later Kormendy & 
Ho (2013) have argued against such conclusion, pointing out that the 
Mgg — My relation holds only for classical bulges and is not followed 
by spirals hosting pseudo-bulges or by bulgeless discs, which would 
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Figure 9. Comparison between our fiducial model (black solid curves in all panels) and models obtained using different parameter values (coloured dashed 
curves). Panels-sets (a) to (e) show the effects of varying a single parameter while keeping the others fixed to their fiducial values. The variations A indicated in 
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suppressed (red); iv) BH feedback is made incapable of quenching star formation (blue). 
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Figure 10. Stochastic realisation of our fiducial model at z = 0 including a Gaussian scatter, with a standard deviation of 0.4, in log( er), log( Merit,o/ Mo) 
and a. The color-coding adopted is the same of Fig. 1. Grey circles are used for systems with My < 5 x 10!! Mo or Mgy < 106 Mo, which do not have a 
counterpart in the observed sample. The black-dashed curves show the model without scatter. 
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Figure 11. Residual correlation plot, analogous to that of Fig. 2, for a stochas- 
tic realisation of our fiducial model at z = 0 including a Gaussian scatter of 0.4 
in log;o( er), log yg (Merit,o/ Mo) and a. Only systems with My > 5x 10!! Mo 
and Mpy > 10° Mo are shown. 


suggest that that it comes as a byproduct of the co-evolution be- 
tween spheroids - arguably formed by mergers and therefore strictly 
related to the accretion history of halos - and BHs. In our study we 
have shown that all galaxies participate in the same Mgy — My re- 


lation (left panel in Fig. 1). This seems to be the case also for the 
12 spirals that host pseudo-bulges, although our statistics are heavily 
affected by the limited mass range where these are present within 
our sample (11.5 < logj9(Mp/ Mo) < 12.3). The Mgy — My rela- 
tion, however, appears to change slope at My of a few x10!2 Mo 
(or Mpy of 107 2 105 Mo), which suggests an origin from different 
competing mechanisms. In the theoretical framework of Section 4, 
these competing mechanisms are the stellar feedback at low Mp and 
the AGN feedback at high My. At My x 1012 Mo, feedback from 
stars regulate the density of cold gas near the BH, and since in our 
model the stellar feedback efficiency depends on Mp via eq. (B4), a 
Mpu — My relation follows. At higher My, the BH regulates its own 
growth as its accretion rate is limited by the balance between the 
AGN feedback energy and the cold gas gravitational binding energy. 
Since the gas cooling rates and the binding energy both depend on 
Mh, another Mgy — My relation follows, albeit with a different slope. 
Ultimately, a relation between Mgy, M, and Mp must be expected 
in any galaxy evolution framework where the material used for the 
stellar and BH build-up is provided by the cooling of the halo gas 
reservoir. 


de Nicola et al. (2019) used a sample of 83 galaxies with different 
morphologies and high-quality Mgy measurements to investigate 
the correlations between Mpy and the host galaxy properties. Not 
surprisingly, our data are in perfect agreement with the best-fit Mgy— 
Ce relation of de Nicola et al. (2019) (dotted black line in the bottom- 
right panel of Fig. 4) as the majority of our c» measurements come 
from their study. More noticeably, they also agree well with the 
relation reported by Kormendy & Ho (2013) for classical bulges 
(red-dashed line in the same panel), especially in the high-mass 
regime. Interestingly, galaxies with oe < 150 — 200 km s7! show a 
characteristic deviation from the best-fit Mpgg — ce line that closely 
mimics that observed in the Mgy — Mp plane. This is only partially 
visible in our Fig. 4, but is much more evident in the larger sample of 
de Nicola et al. (2019, see the top panel of their Fig.2). This points 
towards the existence of a tight correlation between Mp and oe. 


Davis et al. (2019b) studied a sample of 44 spiral galaxies with 
dynamical measurements of Mpy, using the line-width of the in- 
tegrated velocity profile (from Hı or Ha data) as a proxy for vga, 
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which is then converted to Mp using the calibration from Katz et al. 
(20192). A similar study was more recently done by Smith et al. 
(2021) using both spatially resolved and unresolved CO observa- 
tions, and including early-type systems too. The Mpy — My relations 
resulting from these two studies are shown in the top-left panel? of 
Fig. 4, and are consistent with our data only in the low-mass regime. 
However, there are differences between our study and those of Davis 
et al. (2019b) and Smith et al. (2021) in terms of both methods and 
goals. As for methods, we have used a more complete galaxy sam- 
ple spanning different morphological types, and relied on separate 
techniques to determine Mp depending on the galaxy morphology: 
spatially resolved rotation curves in late-types, which we carefully 
selected from the literature, and globular cluster dynamics in early- 
types, mostly coming from the work of PF21. This allowed us to 
see features in the data, namely the break in the relation, that were 
not apparent in the samples of other studies. As for goals, rather 
than focusing on characterising the Mgy — Mp relation in order to 
use one quantity as a proxy for the other, we made the attempt to 
encapsulate our results in a theoretical framework, with the purpose 
of understanding the physical origin for the observed trends. 

In a series of recent works, Davis et al. (2018, 2019a,b) have fo- 
cused on characterising the scaling relations between Mpy and the 
host galaxy properties in late-type systems, which until then had re- 
ceived little attention in the literature. The scaling relations resulting 
from their studies, shown as blue-dashed lines in the various panels 
of Fig. 4, are in good agreement with our data when only late-type 
(Hubble T > 0) galaxies are considered, but are in tension with the 
full sample which includes many early-type objects. This discrepancy 
is particularly severe in the Mgy — Mbuige plane (bottom-left panel 
of Fig. 4), where our data indicate a much shallower slope, virtually 
identical to that found by Kormendy & Ho (2013) for classical bulges 
(red-dashed line in the same panel), with respect to the scaling found 
by Davis et al. (20192). Davis et al. (2018) compared their Mpy — Mx 
relation determined in 40 spirals with that derived using 21 early-type 
galaxies from the sample of Savorgnan et al. (2016), finding signif- 
icant differences in the slope and normalisation. The two relations 
reported by Davis et al. (2018) are shown in the top-middle panel of 
Fig. 4, and indeed appear to better describe separately the two galaxy 
types. Our fiducial model, shown as a magenta dot-dashed line in the 
same panel, is ignorant of galaxy morphology and predicts a unique 
relation that passes mostly in between the two sequences of Davis 
et al. (2018). 

Our theoretical model is largely inspired by the work of B17, who 
however focused almost entirely on how BH growth in halos is related 
to the development of the ‘red’ and ‘blue’ galaxy sequences in the 
present-day Universe. B17 clearly showed the peculiar shape of the 
Mpy - My, relation, which we found to be in excellent agreement 
with the data, and suggested the use of equilibrium models as possible 
improvements for their approach. We have followed their suggestion 
here so that f, can be predicted from within the same framework 
that models the evolution of dark matter halos and BHs. 

While we model the cooling of gas reservoir in halos, it must be ex- 
pected that BH fuelling is primarily governed by physical processes 
occurring at sub-kpc scales rather than by large-scale cosmological 
infall. Recently, Hopkins et al. (2021) presented a model in which 
the fraction of gas available for BH fuelling is regulated by feed- 
back from star formation occurring in the central (R « 1 kpc) galaxy 


3 Here we used the vg,;-to- Mf calibration from eq. (A1) and, for the results 
of Smith et al. (2021), assumed vfat ~ Wso/2sin(i), i being the galaxy 
inclination 
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regions. At these radii, the timescales over which stellar feedback 
operates are longer than the dynamical timescales: this prevents the 
onset of an ‘equilibrium’ phase where the thickness of a the gas 
layer smoothly adjusts to the turbulence injected by stellar feedback 
(e.g. Marasco et al. 2015; Bacchini et al. 2020), producing instead 
a galaxy-scale outflow with an efficiency that depends on the cen- 
tral surface density. This model implies that BH masses trace the 
host galaxy properties above a critical surface brightness (typical of 
bulges), and correctly predicts the slope and normalisation of the 
Mpu-o and the Mgy-Mbuige relations. Our model is not tailored 
to reproduce ‘local’ galaxy properties but it includes a conceptually 
similar prescription, inherited from B17: the density of gas near the 
BH scales as the inverse of the mass loading factor £ (eq. B8). Since 
B depends on Mp (eq. B4), galaxies with My < Merit feature a severe 
deficit in their central gas density and their BHs grow slowly. Con- 
versely, stellar feedback has virtually no impact on the BH growth at 
higher Mp, which is regulated instead by the balance between Epy 
and E54! (unlike the model of Hopkins et al. 2021, which does not 
consider AGN feedback). 

We stress that more refined semi-analytical models than the one 
presented here are available in the literature (e.g. Croton et al. 2006; 
Somerville et al. 2008; Guo et al. 2013; Henriques et al. 2015; 
Behroozi et al. 2019). These approaches include several ingredients 
such as mergers, environmental effects, different modes of AGN- 
driven feedback, and are designed to capture the evolution of chemi- 
cal abundances and angular momentum of galaxies in addition to the 
quantities studied here. However, none of these works have focused 
specifically on the Mpy — Mp — fx relation, most likely because of 
the lack of high-quality data available at the time of the writing. We 
expect that the dataset built in this work (Table 1) will be useful for 
theorists to constrain future models of galaxy evolution. 

Cosmological hydrodynamical simulations in the ACDM frame- 
work such as Illustris (Vogelsberger et al. 2014) predict relations be- 
tween My, Mpy, Mp and spiral arm pitch angle in qualitative agree- 
ment with the observations (Mutlu-Pakdil et al. 2018). A more quan- 
titative comparison with between these predictions and the dataset 
built in the present work would provide useful constraints to galaxy 
evolution models. 


5.2 Considerations on galaxy morphology 


It is well known that galaxy colour and morphology are strongly 
correlated: while the galaxy population transits from the *blue cloud* 
of star forming systems to the ‘red sequence’ of quenched ones, it 
also undergoes a structural transformation from a preferentially disc- 
dominated to a more spheroid-dominated morphology (Roberts & 
Haynes 1994; Baldry et al. 2004; Muzzin et al. 2013; Kelvin et al. 
2014). In this Section we explore the morphology distribution of our 
sample and interpret it in the framework of our model. 

We turn our attention to the color-coding used in Fig. 6, which indi- 
cates the morphological Hubble T-type. The morphology distribution 
is strongly bi-modal in Mj, with late- (early-) type systems systemati- 
cally occupying halos with masses below (above) ~ 3x 10? Mo. This 
corresponds to a similarly sudden transition around fẹ of 0.15—0.2. 
The same bimodal behaviour was already noticed by PF21. The tran- 
sition as a function of Mgy appears instead to be more gradual and 
occurs at BH masses between 107 Mo and 3 x 108 Mo, roughly 
where the break in the BH scaling relations appears. 

We already touched upon the different scalings followed by the 
different morphological types in Section 5.1. Considering earlier and 
later types separately in Fig. 6, one may argue that the two categories 
obey different scaling laws, as the separation in the Mgy-My, space 


and, above all, in the f,-Mpy space is visible. However, our theoreti- 
cal model predicts a unique relation in the Mpg- Mp- fx space. While 
information on morphology (or angular momentum) are not specif- 
ically encoded within our model, one may question the existence of 
a single evolutionary sequence valid for all galaxy types. 

Fully addressing the question above goes beyond the purpose of 
this study, but a qualitative argument comes from the fact that, at a 
fixed Mpy, galaxies with higher fx do not only live in lighter halos 
(see Section 3 and Fig. 2) but also have later morphological types 
(rightmost panel of Fig. 6). In our theoretical framework, lighter ha- 
los at z 2 0 are always associated to later seeds which, due to their 
relatively short lifespan, had a lower probability of experiencing ma- 
jor mergers that could induce a morphological transformation (e.g. 
Hopkins et al. 2010; Martin et al. 2018). Hence our model, which re- 
produces the correlations shown in Fig. 2 (Section 4.2), can also qual- 
itatively explain the trends with the morphological type. We stress 
that, within this framework, the relation between galaxy quenching 
and morphological transformation is indirect: quenching is caused by 
BH growth and feedback, which is more significant in more massive 
halos; more massive halos are produced by earlier seeds, which have 
a higher probability of experiencing morphological transformations 
induced by major mergers. In a more realistic scenario, however, it 
is likely that galaxy mergers have an impact on the growth of the 
central BH. For instance, mergers can funnel gas towards the central 
regions of the galaxy, thus promoting BH growth and triggering AGN 
feedback (e.g. Hopkins et al. 2006). Also, variations in the angular 
momentum of a galaxy have an impact on the stability of its disc 
and can possibly affect the evolution of f, as suggested by Romeo 
et al. (2020). In general, the BH accretion rate may increase if the 
overall angular momentum of the galaxy gets reduced (which may 
be the case after a major merger). These considerations suggest that 
the relationship between galaxy structure and quenching are more 
complex than suggested by our simple framework. 


5.3 Limitations of our model 


One of the main conceptual weakness of the model discussed in Sec- 
tion 4 lies in its lack of an explicit, physically motivated prescription 
for star formation processes caused by the absence of a dedicated 
treatment of the ISM. Star formation is instead fully described by 
eq. (B3), which simply states that the material that does not feed the 
BH gets partitioned into stars and recycled gas. In other words, the 
BH is the preferential target for gas accretion, and only the cold gas 
spared by BH feeding - if AGN feedback permits - participates in star 
formation processes. While this prescription simplifies significantly 
the treatment of star formation in our model, it does not capture the 
commonly accepted idea that gas cooling from the circumgalactic 
medium is firstly deposited onto the ISM by different mechanisms 
(e.g. Marasco et al. 2012; Fraternali 2017), and only later can par- 
ticipate to star formation and BH feeding processes depending on 
the complex physics of this medium. Also, since in our model we 
do not follow gas dynamics explicitly, we cannot address the ques- 
tion on how gas is transported from the scales of the circumgalactic 
medium (tens or hundreds of kpc) down to the scales of BH accretion 
discs (~ 107? po). Dedicated analytical and numerical models are re- 
quired explore the physical processes responsible for gas transport to 
different scales (e.g. Hopkins & Quataert 2011). 

The adopted prescription also does not allow us to follow the 
physics of star formation quenching, which is instead implemented 
*ad-hoc' in the model: the time when the BH enters the self-regulated 
phase marks the end of the stellar mass built-up, as any subsequent 
inflow of gas will target the BH alone. This has a number of side 
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effects. First, it produces a net segregation between star-forming and 
quenched systems (Fig. 7), leaving no space for a transitory, ‘green- 
valley’-like regime. Second, as the BH population keeps growing in 
mass while star formation is stopped, the Mp -to- M, ratio increases 
with time in most massive systems (panel d in Fig. 8), which is in 
tension with high-z observations (e.g. Merloni et al. 2010). Third, 
it leads to a SHMR that is slightly too steep in the high-mass end. 
This is only marginally visible from Fig. 6 but becomes much more 
evident when an isothermal equation of state (y = 1 in eq. (B6)) is 
used, as discussed in Section 4.2. A model in which gas accretion 
gets partioned a-priori between star formation and BH feeding, in 
the delicate phase where Egg = E991, may help to alleviate these 
issues and would likely provide a more realistic picture of stellar and 
BH mass assembly in high-mass galaxies. 

Our model assumes a continuous and smooth mass assembly, 
which begs the question if mergers play any role in it. The smooth 
accretion of dark matter and gas that we adopt (eq. B1) can also be 
interpreted as a continuous ‘rain’ of gas-rich minor mergers, but we 
stress that in our framework the bulk of stellar and BH mass is built 
‘in-situ’. While cosmological hydrodynamical simulations suggest 
that most of the stellar mass in halos with Mp > 10? Mo is built 
ex-situ (Pillepich et al. 2018b), which may also explain why Mpuige 
correlates with Mpy better than Misc, a scenario where the BH scal- 
ing relations are preferentially produced by mergers does not seem 
to be favoured by the data (King & Nealon 2021). However, major 
mergers may be relevant to set galaxy morphology, as discussed in 
Section 5.2. 

Finally, we notice that the BH feedback efficiency adopted in our 
fiducial model, ef = 1.0 x 1072, is 5 — 10 times smaller than the typ- 
ical values used in semi-analytical models (e.g. Croton et al. 2006, 
2016) and cosmological hydrodynamical simulations (e.g. Di Matteo 
et al. 2005; Schaye et al. 2015; Pillepich et al. 2018a) to reproduce 
galaxy scaling relations. It is not trivial to compare our feedback pre- 
scription with those adopted in these more complex models, which 
often employ different feedback ‘modes’ (e.g. ‘quasar’-mode and 
*radio'-mode) depending on the BH accretion rate. Also, eg varies 
significantly from one model to another: for instance, in their hy- 
drodynamical simulations, Valentini et al. (2020) implemented the 
coupling between the AGN energy and a realistic multi-phase ISM, 
finding that ep ~ 107? adequately regulates the BH growth in Milky 
Way-like galaxies. A simplification that we adopted here is the ab- 
sence of an AGN duty cycle: our model assumes a continuous feeding 
of the BH, corresponding to a continuous energy output distributed 
over the entire galaxy lifespan, whereas the observed AGN activ- 
ity is intermittent. We plan to take this effect into account in future 
developments of our model. 


6 CONCLUSIONS 


Galaxy formation is regulated by the competition between processes 
that favour the growth of stars and supermassive BHs, such as the 
cooling and gravitational collapse of gas accreted from the intergalac- 
tic medium, and those that quench them, namely negative feedback 
resulting from star formation and BH accretion. The stellar, halo and 
BH masses of present-day galaxies are precious remnants of these 
evolutionary processes and of the key parameters that describe them. 
For these reasons, the characterisation of the relations between Mpy, 
My and M, (or equivalently, as in this study, fx) at any redshift is a 
subject of great interest for both the observational and the theoretical 
astrophysical communities. 

In this work we have focused on the relations between Mgy, Mp 
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and fy atz ~ O using a sample of 55 galaxies with dynamical 
estimates of their BH and halo masses, all having Mgy > 10° Mo 
and My > 5 x 10!!Mo. Unlike previous studies, where the line- 
widths of the integrated velocity profiles from H1, CO or Ha data 
were used as proxies for My (e.g. Sabra et al. 2015; Davis et al. 201 9b; 
Smith et al. 2021), our sample is almost entirely made of galaxies 
for which either spatially resolved rotation curves or globular cluster 
kinematics were previously analysed in the literature, providing a 
more refined estimate of their Mp. Our main results, valid in the 
mass range considered, can be summarised as follows: 


e My and Mpg strongly correlate with each other and anti- 
correlate with f, (Fig. 1). In general, galaxies tend to follow a one- 
dimensional sequence in the Mp-MBH-fx space, rather then dis- 
tributing across a two-dimensional plane (Fig. 3). 

e In our sample, the strength and tightness of the correlation 
between Mpy and Mp is comparable to the one between Mpy and 
Ge, the mean stellar velocity dispersion within the effective radius 
(Fig. 4). 

e Bulge masses correlates with Mpy significantly better than disc 
masses, but not statistically better than the total M, of the host galaxy. 
The limitations of our sample prevent us from assessing whether this 
property is also applicable to a complete, volume-limited sample. 

e The slopes of the Mgu-Mp and the f,-Mpy relations show a 
break at Mp ~ 10? Mo or Mpy ~ 107 — 105 Mo (Fig. 1). 

e Atafixed Mpy, fx correlates negatively with M, and positively 
with the Hubble morphological T-type. That is, at a given BH mass, 
galaxies with a higher global star formation efficiency tend to have 
later morphological types and to occupy less massive halos (Fig. 6 
and 2). There is no significant trend instead between Mpy and f, at 
fixed Mp. 


In Section 5 we developed an equilibrium model in the ACDM 
framework to explain the observed trends and provide insights into 
their physical origin. Our model, largely inspired by that developed 
by B17, assumes that galaxies evolve by smoothly accreting dark 
and baryonic matter at a cosmological rate, while the competition 
between the cooling of the available gas reservoir and negative feed- 
back from star formation and AGN regulates the growth rates of stars 
and BHs. The model is based on five free parameters regulating the 
stellar and AGN feedback, which we manually tuned to match the 
data. In spite of its simplicity, the model reproduces the observed re- 
lations remarkably well, including the break at Mgy ~ 107 21089 Mo 
and the trends found at fixed Mgy and Mp (Fig. 6 and 11). In the 
model, the break originates as the BH population transits from a 
rapidly accreting phase, in which stellar feedback is inefficient and 
the BH feedback energy Egy is still small compared to the gravi- 
tational binding energy of the cooling gas Ecool, to a more gradual 
and self-regulated growth, driven by the continuous balance between 
the Egg and Ego}. This balance produces a slope in the Mpy— My 
relation in excellent agreement with the data at Mp > 10? Mo. The 
correlations at fixed Mpy are instead produced by scatter in the BH 
feedback efficiency: higher efficiencies lead to less massive BHs at 
a given halo mass but do not alter the stellar-to-halo mass relation; 
hence, at fixed Mgy, galaxies with a higher fstar will be those occu- 
pying less massive halos. 

As our model lacks a dedicated treatment of mass accretion pro- 
vided by mergers, it does not allow us to follow separately the forma- 
tion of spheroids and discs, which in turns prevents us from drawing 
conclusions on the physics that regulate the co-evolution between 
spheroids and BHs. More advanced models are required to tackle 
this topic, which we plan to implement in a follow-up study. 
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APPENDIX A: THE CONVERSION FROM Vg Ar TO My 


In Section 2 we make use of rotation curve measurements to estimate 
Mh in our systems. Rather than resorting to theoretical models, we 
prefer to calibrate the My — vga relation empirically using the results 
of Posti et al. (20192), which are based on mass decomposition of 
rotation curves from the SPARC galaxy sample (Lelli et al. 2016). 
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Figure A1. Relation between Mp and vfat for a sample of 125 disc galaxies 
from the SPARC sample (Lelli et al. 2016), as derived from the mass mod- 
elling of their rotation curves from Posti et al. (20192). The shaded region 
shows a linear fit to the data and its rms scatter around it, and is used as a 
calibrator for the vga;-to- My conversion in our sample. 


In Fig. Al we show that this relation holds over a large dynamical 
range, and is given by 


My Vflat 
logio (s = (5.9340.02) + (2.65+0.07) logio | Ec , (Al) 
with an rms scatter of 0.24 dex, which we adopt as the uncertainty 
associated to our Mp estimates via eq. (A1). The final uncertainty 
also incorporates measurement errors in vgar. These, however, are 
often very small or not quoted at all in the literature. 

We note that similar relations for the SPARC sample were recently 
derived by Katz et al. (2019b) using a variety of different assump- 
tions in the rotation curve decomposition procedures. In Katz et al. 
(2019b), the values for the slope (intercept) of eq. (A1) vary in the 
range 2.2—2.9 (5.3 — 6.9), depending on the halo profile adopted and 
on the priors used in the modelling. Our estimates in eq. (A1) are 
well within these ranges. 


APPENDIX B: DETAILED MODEL DESCRIPTION 


In this Appendix we describe in detail our theoretical model of 
galaxy evolution used in Section 4. For consistency with the ob- 
served dataset, we assume that all halo ‘virial’ quantities (such as 
mass, radius and velocity) are defined at (or within) pe. the ra- 
dius where the mean halo density becomes equal to 200 times the 
(redshift-dependant) critical density of the Universe, and simply use 


the suffix *h' to refer to halo properties. 


Halo growth 


We start by assuming a cosmological growth for dark matter halos 
as parametrised by Correa et al. (2015)*: 


Mn 


© 


h 3/2 
"no [70.24 + 0.75(1 * z)] Az 


^ Note that this parametrisation is valid for Mp = M53. as we assume here. 


(B1) 


where ho; = Ho/(70kms~! Mpc™!), Ho is the Hubble constant, 
and Az = (Qn(1 +z)? + Q4)1/3 (so that the Hubble parameter H(z) 
is given by HoA). 

Eq. (B1), like the other equations in this Appendix, is solved nu- 
merically from f = tseeq to t = 13.8 Gyr, corresponding to the cur- 
rent age of the Universe, using a constant time-step ôt of 10 Myr 
and assuming Mp(t = ted) = 101 Mo. By varying teeeq from 
0.3 Gyr to 4 Gyr we can sample halo masses at the present epoch 
from 10105 Mo (for younger seeds) to 1014 Mg (for older seeds). 
We note that varying the seed halo mass is equivalent of varying the 
seed time, and simply allows us to modify the range of Mp sampled 
at z = 0 with no consequences on our results. 

Halos accrete baryons as well as dark matter. At any given time 
t, the total baryonic content of our system is always fp Mp(t), fo 
being the universal baryonic fraction. This material is divided into 
gas, stars, and BH as described below. 


Gas cooling 


We assume that, at a any given time f, the total gas reservoir of our 
model galaxy is given by Mgas(t) = fiMn(t) - My(t) - Mpu(t), 
M,(t) and Mpgg(f) being the stellar and BH mass of the galaxy. 
In reality, this gaseous reservoir is a complex, multi-phase structure 
build by a mixture of heterogeneous processes, including cosmolog- 
ical accretion of pristine gas, outflows produced by stellar and BH 
feedback. Our model does not aim to capture the complex physics of 
this medium. We make instead the simplifying assumption that the 
entire reservoir is in a single phase defined by the virial temperature 
Th of the dark matter halo. Therefore, in what follows we make no 
distinction between interstellar and circumgalactic media, and sim- 
ply refer to the whole gas reservoir as the ‘gas’ component of the 
model. 

Clearly, only a small fraction of Mgas can be used to form stars 
and to feed the BH, and we assume that this fraction is determined 
by the cooling rate of the gas reservoir. Specifically, following White 
& Frenk (1991), at any time t we define the halo cooling radius rego] 
as the radius within a halo where the cooling time is equal to t. Using 
the formulation from Chen et al. (2020), we have 

1/2 
Maas A (Tt. Z) | (B2) 


Teool = MIN 4 rh, ete 
4nump(5kBThvh) 


where ry, Th and vy are the halo virial radius, temperature and circular 
velocity measured at rp, respectively, A(T}, Zgas) is the cooling func- 
tion (Zgas being the gas metallicity), kg is the Boltzmann constant 
and mp is the proton mass. In this work we use the collisional ioni- 
sation equilibrium cooling functions of Sutherland & Dopita (1993) 
at metallicity Zgas = 0, representative for a pristine gas. In Section 
4.2 we explore the effect of using a higher metallicity. 

The ratio reoo1/7p is a decreasing function of both cosmic time 
and My. At high z, where halo masses are small, densities are high 
and gas can efficiently cool, rgoo)/7p = 1 in all systems. As time 
progresses, halos grow and cooling remains efficient only within a 
fraction of the virial radius, but in the lowest mass regime (Mp(z = 
0) € 2x 10!! Mo) the ratio remains 1 even at z=0. This is related 
to the dichotomy between the so-called ‘cold’ and ‘hot’ mode gas 
accretion (e.g. Kere§ et al. 2005; Dekel & Birnboim 2006), with the 
former (latter) working preferentially at high (low) z and in low- 
mass (high-mass) systems. In present-day Milky-Way-like galaxies, 
cool /Th decreases below 1 at z ~ 4. 
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Assuming NFW profiles and the time-dependent M,- 
concentration relation from Dutton & Macció (2014), we can com- 
pute fooo|, the mass fraction enclosed within roo], and the total 
mass of the cooling gas Mgoo1 as feoo!Mgas- The time-derivative of 
this quantity, Mcool, gives the rate that becomes gas eligible for star 
formation and BH feeding at any given time. 

We implicitly assume that the dominant timescale of the gas ac- 
cretion process is the gas cooling time (e.g. Binney et al. 2009), and 
that the newly cooled material can be used immediately by the galaxy 
with no delay time, that is tayn « fcogl. 


Stellar growth and feedback 


We assume a simple equilibrium model where the the rate at which 
the gas cools is balanced by the rates at which stars form, the BH 
grows and gas is returned to the initial reservoir because of stellar 
mass losses and feedback. Such a model can be described by 


Moo! = Mpu + (1 -R + B)M. (B3) 


where & is the recycled gas fraction due to stellar mass losses and £ 
is the mass-loading factor of the outflows driven by stellar feedback. 
Here we set R = 0.3 (e.g. Fraternali & Tomassetti 2012), but similar 
results can be obtained using R = 0.5 (Bruzual & Charlot 2003) by 
slightly adjusting the model parameters. The parameter f, instead, 
requires a more in-depth discussion. 

Following an approach that is very popular in semi-analytical mod- 
els, and is supported by arguments based on energy- and momentum- 
driven winds, we assume that £ is a decreasing function of Mj (e.g. 
Somerville & Davé 2015, and references therein). This choice leads 
to a scenario where stellar feedback is very efficient in the low-mass 
regime, leading to low star formation efficiencies, but becomes less 
effective at larger masses. We parametrise the mass-loading factor as 


7 Mp =g 
j z (x) BS 


with a > 0. In eq.(B4), Merit is a time-dependent critical halo mass 
below which stellar feedback is more effective. As in B17 and Dayal 
et al. (2019), we assume that Merit scales with redshift as 


—3/8 
Merit = Merito Az! (B5) 


where Merit,0 sets the critical mass at z = 0. We notice that the redshift 
dependence in eq. (B5) is very weak: me =1 at z=0, ~ 0.76 at 
z=2 and ~ 0.59 at z = 5. In B17, æ is fixed to 8/9 and Merit,o to 
10? Mo, while Dayal et al. (2019) adopt Merito 101125 Mo. In 
this work, we prefer to treat both œ and Merit,o as free parameters. 

Eq. (B3) can be used to compute the stellar mass growth in the 
time interval [f, t + ôt] if the BH mass growth is known. We de- 
scribe how the latter is determined below. While a starting seed 
Myx seed = Mx(t = fseeq) is required to integrate eq. (B3), in prac- 
tice we find that systems quickly lose memory of their initial seed 
and converge towards a stable solution for any reasonable choice 
of My seed- Mx,seed has some impact on systems with present-day 
My, < 10195 Mo but, in the range of masses studied here, can be 
considered only a technical necessity due to our implementation and 
has no particular physical meaning. For simplicity, we set M, seed to 
10? x MBH, seed- 

In our model, feedback-driven outflows (via £9) and stellar mass 
losses (via R) immediately re-join the gas reservoir (so that the 
total baryonic mass within the virial radius remains fp Mp) and can 
participate in future stellar and BH growth. 
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BH growth and feedback 


We assume Bondi-like accretion (Bondi 1952), where the BH grows 
at a rate 


(B6) 


where G is the gravitational constant, ppgg is the density of the gas 
near the BH (ideally measured at the Bondi radius) and c, is the 
gas effective sound speed. Eq. (B6) represents a specific expression 
to parametrise the BH accretion as a function of BH mass and the 
thermo-dynamic properties of the gas, and has been largely adopted 
in the literature using different normalisation factors (e.g. Di Matteo 
et al. 2005; Croton et al. 2006; Bower et al. 2017). We stress that 
different prescriptions are also possible (e.g. Hopkins & Quataert 
2011), and that the mode by which BHs grow in the early Universe 
can be significantly different (e.g. Kroupa et al. 2020). 

By adopting a polytropic equation of state Pay c Plus Ppg being 
the gas pressure, eq. (B6) can be written as 


im F 3/2 
T ce 2 p Peos 2. (5-3y)/2 
MBH 4nG y kp Teos | MH PBH (B7) 


where u = 0.6 is the gas mean atomic weight, mp is the pro- 
ton mass, kg is the Boltzmann constant, and peos and Teos define 
the normalisation of the equation of state. Here, for simplicity, we 
adopt Teos = 8000 K and Peos/HMp = 0.375 cm~, meaning that, at 
the temperatures typical for the warm ISM, the pressure Peos/kp 
is 3000 Kcm73, similar to that determined for the Galactic ISM 
(Wolfire et al. 2003). However, we stress that any choice for Teos 
Or peos is degenerate with the normalisation of pp defined below 
in eq. (B8), thus the exact values of these quantities cannot be con- 
strained by our analysis. For consistency with B17 we adopt y = 4/3, 


so that Mgy « Miu Pug . Thus the BH growth will strongly depend 
on the BH seed mass, which is another free parameter of our model, 
and weakly on the gas density, which we model as described below. 
The isothermal case (y = 1) is briefly discussed in Section 4.2. 

As in B17, we make the assumption that the density of the gas near 
the BH scales proportionally to the typical density of the ISM, and is 
inversely proportional to the mass loading factor £ as feedback from 
star formation is supposed to remove gas from the innermost regions 
of the galaxy. As discussed, in this study we do not model the ISM 
explicitly, but we can make an educated guess on how its density 
scales with the main ingredients of our model. By assuming that the 
ISM stratifies in a disc with a radially constant thickness, we write its 
density as cc MīsM/ [NE Simple treatments for these quantities are 
Mjsm « Mhp and, in a scenario where the halo angular momentum 
sets the size of the disc (Fall & Efstathiou 1980; Mo et al. 1998), 
Rism & rh & Mi a7! (see Posti et al. 2019b, for a more detailed 
treatment of disc sizes). With these assumptions, ppy can be finally 
written as 


l 
PBH ^ PBH,0 lw BN (B8) 
9 A102 Mg 


where the normalisation ppp, o is a free parameter. 

Our expression for pgp follows the prescription of B17 and is 
meant to capture plausible variations in the central gas density, but 
we stress that different prescriptions can be adopted. For instance, we 
experimented using MisM « Moo, finding results perfectly compat- 
ible with our current implementation. We also tried a more straight- 
forward approach with pgg œ Mooo) finding similar results, except 
for a somewhat larger slope in the SHMR in the high-mass range. All 
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considered, our results do not depend significantly on our prescrip- 
tion for ppp. 

As M, grows with time, feedback from star formation becomes 
increasingly less efficient (eq. B4) and Mgy would rapidly diverge 
if there were no mechanisms that prevented inflow, such as AGN 
feedback. In our model, we implement a ‘preventative’ AGN feedback 
recipe (e.g. Zinger et al. 2020), where the energy released by the BH 
counteracts gas cooling and accretion onto the galaxy. 

A BH radiates away a fraction e; of its accreted rest-mass energy 
providing a luminosity LAGN = €Mpg c?,c being the speed of light. 
Since only a fraction e¢ of such luminosity couples with the gas, the 
AGN feedback energy rate Egy will be 


Égg = €t Laon = & et Mpy c? . (B9) 


As commonly done in the literature, we set e; —0.1 (e.g. Soltan 1982; 
Marconi et al. 2004; Davis & Laor 2011) and treat er, the so-called 
‘BH feedback efficiency’, as a free parameter of our model. 

The BH heats the surrounding gas at a rate given by eq. (B9). 
In principle, a detailed calculation of the balance between heating 
and cooling rates would be required to model self-consistently the 
physics of gas accretion onto the galaxy. However, this is not a trivial 
task, given that we ignore both the scale over which the AGN energy 
is transferred to the gas and the way by which baryons redistribute 
within the halo after the heating and cooling processes. We therefore 
adopt a simplified approach, similar to those often adopted in the 
literature (e.g. Bower et al. 2017; Chen et al. 2020), that aims at 
capturing the effects of this balance integrated over the entire galaxy 
lifespan. We assume that the energy output specified by eq. (B9) is 
accumulated during the lifespan of the BH with no significant energy 
loss. At each time-step, the BH energy output accumulated up to that 
epoch, Epy(t), is compared to the gravitational binding energy of 
the gas within the cooling radius of the system, given by 

cool 

Ecog] = 41G f Pgas(r) Mp(< r)rdr, (B10) 
where pgas (r) is the gas density profile, which we assume to be NFW- 
like, and My(« r) is the enclosed halo mass. We impose Epy(t) < 
Eoo (t), so that the implemented mass accretion onto the BH is given 
by the minimum between eq. (B6) and the value required for Egy(t) 
to be equal to Eggo) (t). The physical justification for this choice is 
that, at any given time, the shell of cooling circumgalactic gas that 
produces an accretion rate equal to Mcoo] is located - by definition - at 
r = Fcool; thus the BH needs to heat all baryons within this radius in 
order to prevent the cooling material from accreting onto the galaxy. 
Since Egy is an integrated quantity, at early times Egy « Eco. 
and the BH population can grow unimpeded following eq. (B6) and 
(B8). Growth rates decrease dramatically when BHs enter the self- 
regulated, binding energy-limited accretion phase. Examples of BH 
mass built-up are given in Section 4.1. 

A key choice that we make is that, at times when the BH is in the 
self-regulated phase, we also stop the stellar mass growth. In other 
words, AGN-driven feedback in our model has two main effects: 
it quenches star formation and regulates BH accretion. In principle 
this phase is not permanent, given that at subsequent times the halo 
accretes new material and the gas cooling proceeds, effectively in- 
creasing Eo}. In practice though, as we discuss in the main text, we 
find that once a BH enters the self-regulated phase it never leaves 
it, and star formation remains suppressed for the rest of the galaxy 
lifetime. 
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